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ABSTRACT 


In this study the effective elastic modulus of cellulose nanocrystal (CNC) 
network is evaluated using multiscale method and micromechanical analysis. 
For this purpose, the elastic modulus of CNC-water phases are randomly 
assigned to a two-dimensional (2D) checkerboard structure and the elastic 
response is evaluated. In addition, the effect of having a different number of 
phases (CNC, water and interface) is evaluated by assigning a discreet and 
continuous distribution of elastic modulus to checkerboard structure. When 
the number of phases increases dramatically, the distribution of phases is 
continuous and is defined with Weibull distribution. The results show that for 
two-phase materials (CNC and water) when the microstructure has organized 
pattern the rule of the mixture and numerical model provide the same 
effective modulus, however when the microstructure is completely random, 
the self-consistent micromechanical model should be used. Also, this study 
suggests 50% volume fraction as the percolation threshold for the CNC 
network with 10 GPa effective elastic modulus. The results from percolated 
multiphase network reveal that for microstructures with 4 phases and above, 
the percolated network converge to 35 GPa. 
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INTRODUCTION 

Cellulose is the most abundant polymer on the planet [1,2]. 
In the past decade new particles of cellulose called 
crystalline nanocellulose (CNC) with outstanding mechanical 
properties received tremendous attention for packaging and 
reinforcement agent [1-3]. Mechanical properties of CNC 
solutions depends on many factors such as the dispersion 
and orientation of CNCs in water which can be solved based 
on percolation theories. When the material reaches the 
percolation threshold (the minimum volume fraction of the 
particles that percolation occurs), the microstructure 
significantly changes and as a result of the macroscopic 
properties significantly change [4]. Percolated CNC network 
could be assumed to be a heterogeneous material due to 
different mechanical properties caused by their density and 
orientation in the network. Mechanical properties of 
heterogeneous materials based on their distribution type can 
be categorized into two general group, (1) materials with 
lognormal distribution of their mechanical properties found 
at some length scale of bone, concrete and many other 
materials and (2) materials with multi-modal distribution 
found in collagen, rocks and polycrystalline materials [5,6]. 
Having a different type of distribution for their mechanical 
properties makes the prediction of their overall mechanical 
properties challenging. One of the models that have been 
employed to resemble a simple microstructure for 
heterogeneous materials is called a checkerboard material 
[5,7]. This type of microstructure is simply formed by 
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assigning different mechanical properties to each partition of 
a 2D partitioned square. The benefit of using such structure 
is that, due to having many partitions in the geometry, by 
assigning different mechanical properties one can form 
multi-phase materials [5]. For example, by randomly 
assigning two elastic modulus one can form randomly 
distributed two-phased materials. Majority of the research 
on the checkerboard structure are related to effective 
conductivity [7-9] and very few on elastic modulus [5]. 
Dimas et al. (2015), studied how the effective stiffness of 2D 
checkerboard affected by the size of tiles and modulus 
contrast between phases and obtained analytical 
approximations based on the parameters of their study on 
percolation of two-phases materials [5]. In addition to 
numerical models, micromechanical models can provide the 
effective mechanical properties of a complex system such as 
two-phase randomly distributed or multi-phased materials. 
There many micromechanical models developed for specific 
composite materials and may be suitable for a certain case 
and volume fraction or type of distribution, while they may 
fail for other cases [9]. In this work, we use numerical 
checkboard model and two micrometrical model, (1) a 
simple micromechanical model obtained by a mixture of 
Voigt and Ruess model and (2) self-consistent method, to 
characterize the elastic modulus of a percolated CNC 
network based on randomly two-phase and multi-phase 
materials. 
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MATERIALS and METHODS 

We use a multiscale framework to obtain the elastic modulus 
of a percolated CNC network. First, molecular dynamics 
(MD) simulation is used to find the elastic modulus of a CNC 
immersed in water. MD simulations as of the major 
characterizing tools for measuring mechanical, thermal and 
interfacial properties of materials have been used by many 
authors before for CNC and other materials [10-12]. Then, 
we use the finite element method (FEM) and checkboard 
structure to model a percolated system of CNC and water. 
Using multiscale framework helps to reduce the 
computational cost and bridge the scale from nano to macro 
scale as has been done previously by many authors. 
Upscaling the mechanical and interfacial properties from 
Nano to macroscale has been done by many researchers 
before for CNC [13] and carbon nanotube [14]. For example, 
Shishehbor et al. (2018), used a novel atomistic to 
continuum framework for evaluating the role of bonded and 
non-bonded interactions on the elastic properties of CNCs 
and showed covalent bond in the main contributor 
interaction in both bending and tensile loading [13]. We have 
to mention that, in addition to the multiscale framework, 
coarse-grained (CG) molecular dynamics have been used 
extensively to bridge the scale from nano to macro [15-18]. 
For CNC, there are few promising CG models that can bridge 
the scale and perhaps can be used for percolation of CNCs. 
For example, the new CG model developed recently by 
Shishehbor and Zavattieri, (2019) seems to be very 
promising for CNC-water and CNC-composite materials 
[17,19,20]. However, the interfacial energy for CNC and 
water needs to optimized for the model before usage. 

I. CNC-water elastic modulus 

We use MD simulation to find the effective elastic modulus of 
a CNC in a box of water. A box of water with dimensions of 
4.0X4.0X13.5 nm that contains a CNC with dimensions of 
2.5X2.5X12 nm in the middle (as shown in Fig. 1) is created 
using VMD package [21]. The simulation is performed at 
300k temperature, 1 atm pressure using NPT ensemble and 
CHARMM forcefield [22]. All the simulations are performed 
using LAMMPS package [23]. 
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Figurel. Molecular structure of CNC immersed in water. 
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Figure 2. The stress-strain curve at different strain rate 
for CNC-water. (a) the stress-strain curve in the z 
(covalent bond direction) and (b) stress-strain curve in 
the x-direction (van der Waals). 

The stress-strain curves obtained for CNC-water is shown in 
Fig. 2a and Fig.2b in covalent bond (axial) and van der Waals 
(vdW) direction respectively. The result shows that the 
elastic modulus is almost independent of the strain rate and 
equal to 100 GPa and lGPa for covalent bond and vdW 
directions respectively. The elastic modulus obtained here is 
consistent with many MD results obtained previously by 
other researchers [13,17,24], except that due to immersion 
in water, the effective elastic modulus is 40-60% lower than 
170-200 GPa reported values [17]. 


II. Finite element and micromechanical models 

The microstructure of the problem is a 2D square of length 
10, which is meshed with the tile length of 1.0, 0.5 and 0.25 
and provides 10x10,20x20 and 40x40 squares. For the FEM 
analysis, the nodes on left and bottoms are fixed in the x, y- 
direction, the displacement is applied on the right nodes, and 
the nodes on top are free as shown in Figure 3a. The effective 
elastic modulus of the micro structure from FEM analysis is 
obtained by dividing the average stresses by average strains 
calculated in the gauss points of elements. 


^ Ee ff^FEM 
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Where (E e f /)fem is the effective elastic modulus, T av0 is 
the average stress and £ ave is the average strain. For the 

micromechanical model, a two-step averaging method is 
performed. First, the average elastic modulus is obtained in 
each column, - E n , with Voigt model due to the fact that all 

the elements (tiles) in the same column have the same strain. 
Then, the total effective modulus is obtained from Ruess's 
model for all averaged columns have the same stress. 

1 

E\ = — [E\ + E2 + .... E n ] 

" ( 2 ) 

^fEmicro^ 1 = \r[E\ l + E 2 ~ 1 + ...En~\ 

(3) 


Where (E 0 f /) m j wd is the effective elastic modulus from the 

micromechanical model, n is the number of tiles in each 
column and/or row and ff 1 - E n are the effective elastic 

modulus in each column from simple averaging (Voigt 
model)? For two-phase materials, the elastic modulus of stiff 
and soft tiles are i?i = IDO and E 2 = 1 (from MD analysis) 
with equal poison ratio of ^ = $2 =0.3. Here, first, some 
microstructures with specific patterns are studied and after 
verification, the two-phase random microstructure with 
different volume fraction is evaluated. In order to study the 
size effect on the effective elastic modulus, random 
microstructure with different size for tiles (element) is also 
analyzed. 
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Figure3. The microstructure, FEM and micromechanical model, (a) 10x10 mesh and the FEM setups. The displacement is 
applied to the right boundary while the nodes on the left and bottom are fixed in x and y direction respectively, (b) 
Checkerboard microstructure with 50% stiff and 50% soft materials. The average elastic modulus of each column is 
obtained from the Voigt model and denoted as E±- E^ (c) Random checkerboard with 50% stiff and 50% soft matrix. The 
average elastic modulus of each column is obtained from the Voigt model and denoted as E n . 


For multi-phase materials, random distribution of elastic modulus from E^ to Eg is assigned to tiles while the number of phases 
changes from 1 to 10. In addition, a Weibull distribution of elastic modulus with different values for the Weibull slope 
parameter is studied to evaluate the effect of distribution on the elastic modulus. 

RESULTS AND DISCUSSION 

III. Random vs organized microstructure 

To test the validation of the FEM and micromechanical model, different case studies with different micro structure and volume 
fraction is tested. Figure 4. shows the different case studies from very simple microstructure to a completely random case 
where the elastic modulus for black and yellow tile is E ± = 100 and E z = 1 respectively. The effective modulus from FEM and 
micromechanical models are shown below each case study in Figure. 4 denoted as FEM and micro respectively. It can be 
observed that for the cases with certain patterns, the FEM and micromechanical models provide the same results and the slight 
difference is due to the presence of poison effect in FEM and ignoring that in the micromechanical model. Even for a 
checkerboard with 50% stiff tile and 50% soft tile, the micromechanical results is very close to FEM with 15% error. However, 
for a random case study shown in Figure 4.j, the micromechanical model value is 47.5 as opposed to FEM result which is 17.5 
and the error is more than 170%. 
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Figure4. Case studies for different microstructure with different volume fraction and patterns. In all cases, the prediction 
of the micromechanical model is equivalent to FEM results except the random pattern. 

IV. Percolated CNC network 

To better evaluate the random microstructure, different case studies with different volume fraction and randomness are 
studied. Figure 5.a displays the variation of effective elastic modulus with volume fraction for 10x10 and 40x40 squares. Each 
case study was executed 3 times to produce a different random structure. The results indicate that the variation of effective 
elastic modulus with randomness for smaller microstructure, 10x10, is more than 40x40 microstructure with almost the same 
average curve. In addition, comparing the micro and FEM results, it can be realized that the simple averaging micromechanical 
model fails to provide FEM results especially for medium volume fraction (40%-60%). Figure 5b. compares the effective elastic 
modulus from FEM, averaging micromechanical model and self-consistent approach. The results indicate that the self- 
consistent approach provides very accurate results for all range of volume fraction for two-phase random checkerboard 
microstructure. Also, the result shows a sa udden change in effective elastic modulus for 50% volume fraction and therefore, it 
seems that 50% volume fraction with associated 10 GPa effective elastic modulus is the percolation threshold for CNC network. 
Indeed, many studies have estimated the elastic modulus of CNC nanofilm to be around 5-20 GPa depending on the alignment 
and therefore the result shows that this value is in a reasonable range [25]. 



Figure5. Variation of effective elastic modulus with volume fraction, (a) comparing the effect of size of the microstructure 
and randomness on the results, (b) Comparing the effective elastic modulus from FEM, averaging micromechanical model 

and self-consistent approach. 

Two-phase vs multiphase percolated network 

In the previous section, the percolated CNC network was modeled as a two-phase material with E=100 GPa and E=1 GPA. 
However, due to the orientation of CNC particles, each checkboard could take elastic modulus between 1-100 GPa. For multi¬ 
phase materials, the elastic modulus in a range of E1-E2 is randomly assigned to tiles and by controlling the number of phases 
involved in the microstructure, by controlling how many different elastic moduli is provided, the variation of the effective 
elastic modulus with respect to the number of phases can be obtained as shown in Figure 6. When the number of phases 
increases sufficiently, the discrete distribution could be approximated by a continuous one as shown in Figure 6.c. 
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Figure6. Discrete distribution of elastic modulus for having multi-phase materials, (a) Two-phase materials, (b) Three 
phase-materials and (c) Multi-phase materials which can be approximated with continuous distribution. 


Figure.7a displays the variation of effective elastic modulus with a number of phases and its comparison with the 
micromechanical model. It is interesting to see that by increasing the number of phases, the FEM results have an asymptotic 
value. In other words, the effective elastic modulus for having 10 different phases and 5 different phases are almost the same. 
Note that each case study has been performed 4 times to also see the variation of curves with randomness. In addition to 
having a discrete number of phases, having a continuous distribution of elastic modulus provides continuous variation of a 
number of phases. Here, Weibull distribution is used to provide such a variation to the elastic modulus. The equation for the 
probability distribution function of Weibull is: 


/(*) = .£ (£)<*-» e 
XX 


-(t) p 


(4) 


Where (3 is the scale parameter and X is the mean value. The effect of Weibull scale parameter, (3, which controls the diversity of 
the number of phases is also studied as shown in Figure 7b. For (3=1 the FEM and micromechanical model show very different 
results and although the micromechanical value is very close to the average elastic modulus in the Weibull, a=E=50, the FEM 
results are around 30. Comparing the discrete and continuous distribution number of phases shows that when the variation of 
elastic modulus is large, (3=1, the continues distribution is the asymptotic value of discrete number of phases, here around 35, 
but as the (3 parameter increases, even though there are many phases involved in the microstructure due to having continuous 
distribution, the FEM and micro results are close to the average value (around 45). 
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Figure7. The discrete and continuous distribution of a number of phases in the microstructure, (a) The variation of 
effective elastic modulus with respect to a number of phases, (b) The effective elastic modulus for the different scale 

parameter of Weibull distribution. 


CONCLUSION 

In this paper the effective elastic modulus of two-phase and 
multi-phase percolated CNC network was estimated using 
FEM and two different micromechanical models (a mixed 
Voigt-Ruess micromechanical and self-consistent model). 
The microstructure is produced by randomly assigning 
elastic modulus to a 2D checkerboard structure. The results 
show that for two-phase materials when the microstructure 
has a certain pattern the micromechanical model and FEM 


provide the same effective modulus, however when the 
microstructure is completely random, the micromechanical 
model fails to produce FEM results, although using self- 
consistent model shows a great match with FEM results for 
different volume fraction. For random multiphase materials, 
the effect of a number of phases on the effective elastic 
modulus was studied by discreetly controlling the number of 
phases and randomly assigning them to checkerboard tiles. 
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In addition, the effect of a continuous distribution of elastic 
modulus, continuous phases, were studied with Weibull 
distribution of elastic modulus. The results show that by 
increasing the number of phases, the effective modulus for 
discrete distribution has an asymptotic value which is equal 
to Weibull scale parameter (3=1. By increasing the (3, the 
effective elastic modulus gets closer to the single-phase with 
average elastic modulus in the Weibull distribution. In other 
words, having a continuous distribution only provides two 
sides of the discrete distribution, one for two-phase 
materials and one for many phases which provides 
asymptotic effective modulus. Finally, this study suggests 
50% volume fraction as the percolation threshold for CNC 
network with 10 GPa effective elastic modulus. The results 
from percolated multiphase network reveal that for 
microstructures with 4 phases and above, the percolated 
network converge to 35 GPa effective elastic modulus. 
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